home *** CD-ROM | disk | FTP | other *** search
/ APDL Eductation Resources / APDL Eductation Resources.iso / programs / electronic / rlab / TestMatrix / kahan_r < prev    next >
Encoding:
Text File  |  1994-12-20  |  2.0 KB  |  64 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Kahan matrix - upper trapezoidal.
  4.  
  5. // Syntax:    K = kahan ( N , THETA )
  6. //        K = kahan ( N )
  7. //        K = kahan ( N, THETA , PERT )
  8.  
  9. // Description:
  10.  
  11. //    K  is an upper trapezoidal matrix that has some interesting
  12. //    properties regarding estimation of condition and rank. The
  13. //    matrix is N-by-N unless N is a 2-vector, in which case it is
  14. //    N[1]-by-N[2]. The parameter THETA defaults to 1.2. The useful
  15. //    range of THETA is 0 < THETA < PI. 
  16. //
  17. //    To ensure that the QR factorization with column pivoting does
  18. //    not interchange columns in the presence of rounding errors,
  19. //    the diagonal is perturbed by PERT*EPS*diag( [N:1:-1] ).  The
  20. //    default is PERT = 25, which ensures no interchanges for
  21. //    KAHAN(N) up to at least N = 90 in IEEE arithmetic. 
  22. //    KAHAN(N, THETA, PERT) uses the given value of PERT. 
  23.  
  24. //    The inverse of KAHAN(N, THETA) is known explicitly: see Higham
  25. //    (1987, p. 588), for example. 
  26. //    The diagonal perturbation was suggested by Christian Bischof.
  27. //
  28. //      References:
  29. //       W. Kahan, Numerical linear algebra, Canadian Math. Bulletin,
  30. //          9 (1966), pp. 757-801.
  31. //       N.J. Higham, A survey of condition number estimation for
  32. //          triangular matrices, SIAM Review, 29 (1987), pp. 575-596.
  33.  
  34. //    This file is a translation of kahan.m from version 2.0 of
  35. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  36. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  37.  
  38. //-------------------------------------------------------------------//
  39.  
  40. kahan = function ( n , theta , pert )
  41. {
  42.   local (n, theta, pert)
  43.   global (eps)
  44.  
  45.   if (!exist (pert)) { pert = 25; }
  46.   if (!exist (theta)) {theta = 1.2; }
  47.  
  48.   r = n[1];    // Parameter n specifies dimension: r-by-n.
  49.   n = n[max(size(n))];
  50.  
  51.   s = sin(theta);
  52.   c = cos(theta);
  53.  
  54.   U = eye(n,n) - c*triu(ones(n,n), 1);
  55.   U = diag(s.^[0:n-1])*U + pert*eps*diag( [n:1:-1] );
  56.   if (r > n)
  57.   {
  58.     U[r;n] = 0;        // Extend to an r-by-n matrix.
  59.   else
  60.     U = U[1:r;];    // Reduce to an r-by-n matrix.
  61.   }
  62.   return U;
  63. };
  64.